Learningtower: Comparative Analysis of PISA 2022 and Historical Data

Shabarish Sai and Guan Ru Chen

Department of Econometrics and Business Statistics

2024-10-18

Contributors

  • Shabarish Sai Subramanian

  • Guan Ru Chen

  • Dianne Cook

  • Kevin Y.X. Wang

  • Priya Ravindra Dingorkar

Introduction

The learningtower R package is designed to streamline the analysis of OECD’s Programme for International Student Assessment (PISA) data. This package provides access to datasets from 2000 to 2022, allowing researchers to explore trends in education, student performance, and other contextual factors. It simplifies the process of handling large, complex datasets, making it easier to conduct comparative studies across countries and years. Currently, we are updating the 2022 version of the learningtower package to ensure compatibility with the latest PISA data and functionalities

Collection of Data

PISA data is collected every three years from over 70 countries, targeting 15-year-old students. The assessment measures students’ abilities in reading, mathematics, and science through standardized tests. In addition to the tests, questionnaires are administered to students, teachers, and school principals to gather contextual data on educational environments, socio-economic status, and more. This comprehensive approach helps provide insights into factors that affect student performance across different educational systems worldwide.

PISA Dataset

student_data <- readRDS("../Data/student_2022.rds")
print(colnames(student_data))
 [1] "year"        "country"     "school_id"   "student_id" 
 [5] "mother_educ" "father_educ" "gender"      "computer"   
 [9] "internet"    "math"        "read"        "science"    
[13] "stu_wgt"     "desk"        "room"        "dishwasher" 
[17] "television"  "computer_n"  "laptop_n"    "car"        
[21] "book"        "wealth"      "escs"       

The student dataset includes the following columns: year, country, school_id, student_id, mother_educ, father_educ, gender, computer, internet, math, read, science, stu_wgt, desk, room, dishwasher, television, computer_n, laptop_n, car, book, wealth, escs, and curiosity. These columns provide comprehensive details about the students’ background, academic performance, and access to resources, offering a robust dataset for analysis of educational outcomes and socio-economic factors.

Gender Gap Analysis: Maths

# Your data loading, filtering, and transformation logic as provided:
load(here("data/countrycode.rda"))
readRDS(here("data/student_2022.rds")) -> student_2022

# Load the country names, and join
student_country <- left_join(student_2022, countrycode, by = "country")

student <- student_country %>%
  filter(!is.na(gender), !is.na(math), !is.na(stu_wgt))

# Compute average math scores and gender diff
if (!file.exists("data/math_diff_conf_intervals.rda")) {
  math_diff_df <- student %>%
    group_by(gender, country_name) %>%
    summarise(avg = weighted.mean(math, stu_wgt), .groups = "drop") %>%
    ungroup() %>%
    pivot_wider(country_name, names_from = gender, values_from = avg) %>%
    mutate(diff = female - male, country_name = fct_reorder(country_name, diff))

  # Compute bootstrap samples
  set.seed(2024)
  boot_ests_math <- map_dfr(1:100, ~{
    student %>%
      group_by(country_name, gender) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(math, stu_wgt), .groups = "drop") %>%
      pivot_wider(country_name, names_from = gender, values_from = avg) %>%
      mutate(diff = female - male, country_name = fct_reorder(country_name, diff)) %>%
      mutate(boot_id = .x)
  })

  # Compute bootstrap confidence intervals
  math_diff_conf_intervals <- boot_ests_math %>%
    group_by(country_name) %>%
    summarise(lower = sort(diff)[5],
              upper = sort(diff)[95],
              .groups = "drop") %>%
    left_join(math_diff_df, by = "country_name") %>%
    mutate(country_name = fct_reorder(country_name, diff)) %>%
    mutate(score_class = factor(case_when(
      lower < 0 & upper <= 0 ~ "boys",
      lower < 0 & upper >= 0 ~ "nodiff",
      lower >= 0 & upper > 0 ~ "girls"),
      levels = c("boys", "nodiff", "girls")))

  save(math_diff_conf_intervals, file = here("data/math_diff_conf_intervals.rda"))
} else {
  load(here("data/math_diff_conf_intervals.rda"))
}

math_plot <- ggplot(math_diff_conf_intervals,
                    aes(x = diff, y = fct_reorder(country_name, diff), col = score_class)) +
  geom_point(size = 3) +                                           # Larger points for better visibility
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0.5) +    # Error bars with appropriate width
  geom_vline(xintercept = 0, color = "gray", linetype = "dashed") +# Dashed line at x=0 for reference
  scale_colour_manual("", values = c("boys" = "#3288bd",
                                     "nodiff" = "#969696",
                                     "girls" = "#f46d43")) +        # Custom colors for groups
  labs(y = "Country", x = "Difference (Girls - Boys)", title = "Math Gender Gap Analysis") + # Axis titles
  theme_minimal() +                                                # Clean theme
  theme(
    axis.text.y = element_text(size = 8, hjust = 1),               # Smaller y-axis text, aligned left
    axis.text.x = element_text(size = 10),                         # Normal size for x-axis labels
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"), # Centered title, bold
    plot.margin = margin(20, 20, 20, 150)                          # Add margin for better readability
  ) +
  scale_x_continuous(limits = c(-70, 70),                          # Set x-axis limits
                     breaks = seq(-60, 60, 20),                   # Define breaks
                     labels = abs(seq(-60, 60, 20))) +             # Show absolute values for x-axis
  annotate("text", x = 55, y = 1, label = "Girls", size = 4, color = "red") + # Label for Girls
  annotate("text", x = -55, y = 1, label = "Boys", size = 4, color = "blue")  # Label for Boys

# Print the plot
print(math_plot)

Explanation of Gender Gap: Math Scores

With the gender difference in average maths scores (measured as girls’ scores - boys’ scores) on the x-axis, this graphic displays the gender gap analysis in mathematics across several nations. The y-axis lists the countries, and the lines indicate confidence intervals, and each point displays the average score difference. Grey points indicate no discernible gender difference, red points emphasise nations where girls outperform boys, and blue points indicate nations where boys exceed girls. The graph illustrates the different degrees of gender inequality in maths ability, with boys outperforming girls in many nations and the opposite tendency in a small number.

Gender Gap Analysis: Reading Scores

# Your data loading, filtering, and transformation logic as provided:
load(here("data/countrycode.rda"))
readRDS(here("data/student_2022.rds")) -> student_2022

# Load the country names, and join
student_country <- left_join(student_2022, countrycode, by = "country")

# Subset data and drop missing values for reading scores
student <- student_country %>%
  filter(!is.na(gender), !is.na(read), !is.na(stu_wgt))

# Compute average reading scores and gender diff
if (!file.exists("data/read_diff_conf_intervals.rda")) {
  read_diff_df <- student %>%
    group_by(gender, country_name) %>%
    summarise(avg = weighted.mean(read, stu_wgt), .groups = "drop") %>%
    ungroup() %>%
    pivot_wider(country_name, names_from = gender, values_from = avg) %>%
    mutate(diff = female - male, country_name = fct_reorder(country_name, diff))

  # Compute bootstrap samples
  boot_ests_read <- map_dfr(1:100, ~{
    student %>%
      group_by(country_name, gender) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(read, stu_wgt), .groups = "drop") %>%
      ungroup() %>%
      pivot_wider(country_name, names_from = gender, values_from = avg) %>%
      mutate(diff = female - male, country_name = fct_reorder(country_name, diff)) %>%
      mutate(boot_id = .x)
  })

  # Compute bootstrap confidence intervals
  read_diff_conf_intervals <- boot_ests_read %>%
    group_by(country_name) %>%
    summarise(lower = sort(diff)[5],
              upper = sort(diff)[95],
              .groups = "drop") %>%
    left_join(read_diff_df, by = "country_name") %>%
    mutate(country_name = fct_reorder(country_name, diff)) %>%
    mutate(score_class = factor(case_when(
      lower < 0 & upper <= 0 ~ "boys",
      lower < 0 & upper >= 0 ~ "nodiff",
      lower >= 0 & upper > 0 ~ "girls"),
      levels = c("boys", "nodiff", "girls")))

  save(read_diff_conf_intervals, file = here("data/read_diff_conf_intervals.rda"))
} else {
  load(here("data/read_diff_conf_intervals.rda"))
}

# Plot reading scores with ggplot2
read_plot <- ggplot(read_diff_conf_intervals,
                    aes(diff, country_name, col = score_class)) +
  scale_colour_manual("", values = c("boys" = "#3288bd",
                                     "nodiff" = "#969696",
                                     "girls" = "#f46d43")) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0) +
  geom_vline(xintercept = 0, color = "#969696") +
  labs(y = "", x = "", title = "Reading") +
  theme(legend.position = "none") +
  annotate("text", x = 50, y = 1, label = "Girls") +
  annotate("text", x = -50, y = 1, label = "Boys") +
  scale_x_continuous(limits = c(-70, 70),
                     breaks = seq(-60, 60, 20),
                     labels = abs(seq(-60, 60, 20)))

# Convert the ggplot object to an interactive Plotly plot
read_plotly <- ggplotly(read_plot)

# Display the interactive Plotly plot
read_plotly

Explanation of Gender Gap: Reading Scores

An analysis of the gender gap in reading scores across several nations is shown in this graph. The gender gap in average reading scores is shown by the x-axis, which is computed as (Girls’ scores - Boys’ scores). The lines display the bootstrap confidence intervals, and the y-axis lists the nations. Each point on the y-axis reflects the average gender gap in reading performance. The red dots and lines illustrate that, in the majority of countries, girls perform significantly better than boys in reading, with scores veering towards positive values. The global pattern where girls tend to score higher on reading examinations is highlighted by the vertical zero line, which indicates no difference, and the fact that few countries display boys outperforming girls in reading.

Gender Gap Analysis : Science

# Your data loading, filtering, and transformation logic as provided:
load(here("data/countrycode.rda"))
readRDS(here("data/student_2022.rds")) -> student_2022

# Load the country names, and join
student_country <- left_join(student_2022, countrycode, by = "country")

# Subset data and drop missing values for science scores
student <- student_country %>%
  filter(!is.na(gender), !is.na(science), !is.na(stu_wgt))

# Compute average science scores and gender diff
if (!file.exists("data/sci_diff_conf_intervals.rda")) {
  sci_diff_df <- student %>%
    group_by(gender, country_name) %>%
    summarise(avg = weighted.mean(science, stu_wgt), .groups = "drop") %>%
    ungroup() %>%
    pivot_wider(country_name, names_from = gender, values_from = avg) %>%
    mutate(diff = female - male, country_name = fct_reorder(country_name, diff))

  # Compute bootstrap samples
  boot_ests_sci <- map_dfr(1:100, ~{
    student %>%
      group_by(country_name, gender) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(science, stu_wgt), .groups = "drop") %>%
      ungroup() %>%
      pivot_wider(country_name, names_from = gender, values_from = avg) %>%
      mutate(diff = female - male, country_name = fct_reorder(country_name, diff)) %>%
      mutate(boot_id = .x)
  })

  # Compute bootstrap confidence intervals
  sci_diff_conf_intervals <- boot_ests_sci %>%
    group_by(country_name) %>%
    summarise(lower = sort(diff)[5], upper = sort(diff)[95], .groups = "drop") %>%
    left_join(sci_diff_df, by = "country_name") %>%
    mutate(country_name = fct_reorder(country_name, diff)) %>%
    mutate(score_class = factor(case_when(
      lower < 0 & upper <= 0 ~ "boys",
      lower < 0 & upper >= 0 ~ "nodiff",
      lower >= 0 & upper > 0 ~ "girls"),
      levels = c("boys", "nodiff", "girls")))

  save(sci_diff_conf_intervals, file = here("data/sci_diff_conf_intervals.rda"))
} else {
  load(here("data/sci_diff_conf_intervals.rda"))
}

# Plot science scores with ggplot2
sci_plot <- ggplot(sci_diff_conf_intervals,
                    aes(diff, country_name, col = score_class)) +
  scale_colour_manual("", values = c("boys" = "#3288bd",
                                     "nodiff" = "#969696",
                                     "girls" = "#f46d43")) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0) +
  geom_vline(xintercept = 0, color = "#969696") +
  labs(y = "", x = "", title = "Science") +
  theme(legend.position = "none") +
  annotate("text", x = 50, y = 1, label = "Girls") +
  annotate("text", x = -50, y = 1, label = "Boys") +
  scale_x_continuous(limits = c(-70, 70),
                     breaks = seq(-60, 60, 20),
                     labels = abs(seq(-60, 60, 20)))

# Convert the ggplot object to an interactive Plotly plot
sci_plotly <- ggplotly(sci_plot)

# Display the interactive Plotly plot
sci_plotly

Explanation of Gender Analysis: Science Scores

This graph presents a Gender Gap Analysis in science scores across various countries, showing the difference between girls’ and boys’ average science scores. The x-axis represents the gender difference, calculated as( Girl’s scores - Boy’s Scores), while the y-axis lists the countries. The red points and lines indicate that girls outperform boys in science in several countries, while blue points and lines indicate that boys outperform girls. Grey points and lines represent countries where there is no significant gender difference. The vertical line at zero shows no difference, making it easy to see that in most countries, girls tend to perform better than boys in science, as shown by the positive values on the right side of the chart.

theme_map <- function(...) {
  theme_minimal() +
  theme(
    text = element_text(),
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.border = element_blank(),
    panel.grid = element_blank()
  )
}

region2country = function(region_name){
  country_name = case_when(
    region_name == "Brunei Darussalam" ~ "Brunei",
    region_name == "United Kingdom" ~ "UK",
    region_name %in% c("Macau SAR China", "B-S-J-Z (China)",
                        "Hong Kong SAR China") ~ "China",
    region_name == "Korea" ~ "South Korea",
    region_name == "North Macedonia" ~ "Macedonia",
    region_name == "Baku (Azerbaijan)" ~ "Baku",
    region_name %in% c("Moscow Region (RUS)", "Tatarstan (RUS)",
                        "Russian Federation") ~ "Russia",
    region_name == "Slovak Republic" ~ "Slovakia",
    region_name == "Chinese Taipei" ~ "Taiwan",
    region_name == "United States" ~ "USA",
    TRUE ~ as.character(region_name))
}
math_map_data <- math_diff_conf_intervals  %>%
  dplyr::mutate(country_name = region2country(region_name = country_name)) 

world_map <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  fortify() %>%
  rename(country_name = region)

math_world_data <- full_join(
  x = math_map_data,
  y = world_map,
  by = "country_name") %>% 
  rename(Country = country_name,
         math = diff) %>%
  mutate(math = round(math, digits = 2))
read_map_data <- read_diff_conf_intervals %>%
  dplyr::mutate(country = region2country(region_name = country_name))

world_map <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  fortify() %>%
  rename(country_name = region)

read_world_data <- full_join(
  x = read_map_data,
  y = world_map,
  by = "country_name") %>% 
  rename(Country = country_name,
         Reading = diff) %>%
  mutate(Reading = round(Reading, digits = 2))
sci_map_data <- sci_diff_conf_intervals %>%
  dplyr::mutate(country = region2country(region_name = country_name))

world_map <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  fortify() %>%
  rename(country_name = region)

sci_world_data <- full_join(
  x = sci_map_data,
  y = world_map,
  by = "country_name") %>% 
  rename(Country = country_name,
         Science = diff)  %>%
  mutate(Science = round(Science, digits = 2))
math_dat <- math_world_data %>%
  dplyr::select(Country, math, lat, long, group)

read_dat <- read_world_data %>%
  dplyr::select(Country, Reading, lat, long, group)

sci_dat <- sci_world_data %>%
  dplyr::select(Country, Science, lat, long, group)

math_read_dat <- left_join(math_dat,
                           read_dat,
                           by = c("Country","lat", "long", "group"))

math_read_sci_dat <- left_join(math_read_dat,
                           sci_dat,
                           by = c("Country","lat", "long", "group"))

math_read_sci_dat_wider <- math_read_sci_dat %>%
    pivot_longer(cols = c(2,6,7), names_to = "subjects")

mrs_maps <- ggplot(math_read_sci_dat_wider,
       aes(x = long,
           y = lat,
           group = group)) +
  geom_polygon(aes(fill= value,
                   label = Country)) +
  facet_wrap(~subjects, scales = "free", nrow = 3) +
  theme_map() +
  labs(title = "World Map displaying Gender Gap Scores in Math, Reading and Science")  +
  scale_fill_distiller(palette = "Spectral")
ggplotly(mrs_maps, width=800, height=600)

Explanation for Gender Gap by Countries

The picture displays three global maps that illustrate the gender gap scores for three subjects—math, science, and reading—across various geographical locations. The gender gap value is shown by the colour gradient, where greater values are indicated by darker red and lower values by darker green. Each map shows the gender disparity in schooling in different regions, with notable differences between continents. For example, reading displays more red, indicating a wider gender disparity favouring one gender over the other, but maths exhibits green in many places of the world, indicating fewer gender gaps. Though there are some noticeable regional variations, the scientific map looks comparable to the maths map.

EcoSocio Factors

p1 <- ggplot(data = student_country,
             aes(x = math, y = read)) +
  geom_hex() +
  labs(x = "Math Scores",
       y = "Reading Scores") +
  theme(legend.position="none")

p2 <- ggplot(data = student_country,
             aes(x = math, y = science)) +
  geom_hex() +
  labs(x = "Math Scores",
       y = "Science Scores") +
  theme(legend.position = "none")

p3 <- ggplot(data = student_country,
             aes(x = read, y = science)) +
  geom_hex() +
  labs(x = "Reading Scores",
       y = "Science Scores") +
  theme(legend.position="none")

p1+p2+p3

Breakdown of the Plots

  • The first plot (left) contrasts reading scores (Y-axis) with math scores (X-axis). Students with higher arithmetic scores also typically have higher reading scores, as evidenced by the data points’ highest concentration near the centre of both axes.

  • The second plot (middle) contrasts the science scores (Y-axis) with the math scores (X-axis). A link between math and science results is suggested by the distribution’s high density around the centre of both axes, which indicates that students who do well in math also do well in science.

  • The third plot (right) contrasts the science scores (Y-axis) with the reading scores (X-axis). The hexbin illustrates, like the previous graphs, that kids who score higher on reading assessments typically do better in science, with most data points concentrated in the middle.

Temporal Analysis

# Load student data, and filter to country, cache a copy of the data
# to save downloading every time paper is knitted

# if (!file.exists("data/student_all.rda")) {
#   student_all <- load_student("all")
#   save(student_all, file="data/student_all.rda")
# } else {
#   load("data/student_all.rda")
# }
student_all <- load(here("data/student.rda"))
# Give countries their name, subset to four, and select only variables needed
student_country <- left_join(student,
                             countrycode,
                             by = "country") %>%
  dplyr::filter(country_name %in%
                  c("Australia",
                    "Germany",
                    "Peru",
                    "Qatar",
                    "Belgium",
                    "Brazil",
                    "Denmark",
                    "Greece",
                    "Thailand",
                    "Singapore",
                    "Canada",
                    "Portugal")) %>%
  dplyr::select(year, country_name, math, read, science, stu_wgt) %>%
  na.omit() %>%
  pivot_longer(c(math, read, science), names_to = "subject", values_to = "score")

# Compute the bootstrap confidence intervals, and cache result
if (!file.exists("data/all_bs_cf.rda")) {
  all_bootstrap <- map_dfr(1:100, ~{
    student_country %>%
    group_by(country_name, #year,
             subject) %>%
    #sample_n(size = n(), replace = TRUE) %>%
    mutate(year = sample(year, replace=FALSE)) %>%
    group_by(country_name, year,
             subject) %>%
    dplyr::summarise(
      avg = weighted.mean(score, w = stu_wgt, na.rm = TRUE), .groups = "drop") %>%
    #ungroup() %>%
    mutate(boot_id = .x)
  })

  all_bootstrap_ci <- all_bootstrap %>%
    group_by(country_name, year,
             subject) %>%
    summarise(
      lower = min(avg), # sort(avg)[5],
      upper = max(avg), #sort(avg)[95],
      .groups = "drop")

  # compute original estimate of average and join
  all_avg <- student_country %>%
    group_by(country_name, year, subject) %>%
    summarise(
      avg = weighted.mean(score,
                          w = stu_wgt, na.rm = TRUE),
      .groups = "drop")

  all_bs_cf <- left_join(all_avg,
                      all_bootstrap_ci,
                      by = c("country_name",
                             "year",
                             "subject"))

  save(all_bs_cf, file= here("data/all_bs_cf.rda"))

} else {
  load(here("data/all_bs_cf.rda"))
}
all_bs_cf <- all_bs_cf %>%
  mutate(year = as.numeric(as.character(year)),
         country_name = factor(country_name))
                 # levels = c("Singapore",
                 #          "Australia",
                 #          "New Zealand",
                 #          "Germany",
                 #          "Qatar",
                 #          "Indonesia"))


country_names_highlight <- c("Australia", 
                             "Germany", 
                             "Peru", 
                             "Qatar", 
                             "Belgium", 
                             "Brazil", 
                             "Denmark", 
                             "Greece",
                             "Thailand", 
                             "Singapore", 
                             "Canada", 
                             "Portugal")

math_all_bs_cf_plot <- all_bs_cf %>% 
  dplyr::filter(subject == "math") %>% 
  ggplot(aes(x = year, 
             y = avg)) +
  geom_point(alpha = 0.45) +
  geom_line(aes(group = country_name)) +
  gghighlight::gghighlight(country_name %in% country_names_highlight) +
  labs(
    title = "Maths",
     x = "",
     y = "Score") 


read_all_bs_cf_plot <- all_bs_cf %>% 
  dplyr::filter(subject == "read") %>% 
  ggplot(aes(x = year, 
      y = avg)) +
  geom_point(alpha = 0.45) +
  geom_line(aes(group = country_name)) +
  gghighlight::gghighlight(country_name %in% country_names_highlight) +
  labs(
    title = "Reading",
     x = "",
     y = "Score") 

sci_all_bs_cf_plot <- all_bs_cf %>% 
  dplyr::filter(subject == "science") %>% 
  ggplot(aes(x = year, 
      y = avg)) +
  geom_point(alpha = 0.45) +
  geom_line(aes(group = country_name)) +
  gghighlight::gghighlight(country_name %in% country_names_highlight) +
  labs(
    title = "Science",
     x = "",
     y = "Score")


math_all_bs_cf_plot + read_all_bs_cf_plot + sci_all_bs_cf_plot

Explanation for Temporal Analysis

From 2000 to 2022, the three charts show the temporal trends of math, reading, and science student performance scores in various nations. Labels are used to draw attention to certain countries’ performance trends, and each line shows the average score for that nation.

  • Mathematics: Singapore routinely ranks top, whereas Brazil and Peru have lower scores with some positive trends. Around the 500 score point, nations like Belgium, Australia, and Germany continue to perform comparatively steadily.

  • Reading: Australia, Belgium, and Canada continue to do well, while Singapore once again takes the lead. Thailand, Brazil, and Peru perform worse, though they gradually become better.

  • Science: Australia, Germany, and Belgium retain mid-range ratings, while Singapore and Canada perform at the top. Despite having lower scores, Brazil and Peru have shown some development.

Limitations

Although the Learningtower package makes it easier to access the PISA dataset, it has drawbacks, including less customisation options for more complex analyses, performance problems with huge datasets, and possible incompatibilities with other R versions. Furthermore, for more complicated use scenarios, the documentation might not be adequate. Inherent limitations of the PISA dataset include the fact that it is cross-sectional, which precludes longitudinal tracking, the possibility of sample biases, and the difficulties caused by linguistic and cultural differences that could compromise the comparability of results. Furthermore, the depth of analysis may be constrained by incomplete or missing data, a lack of socioeconomic indicators, and out-of-date background questionnaires. Lastly, the emphasis on standardised test scores may obscure more important educational objectives that the tests do not measure, such creativity and critical thinking.

Thank You